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Abstract 

Recently, it has been proven that evolutionary algorithms produce 
good results for a wide range of combinatorial optimization prob- 
lems. Some of the considered problems are tackled by evolutionary 
algorithms that use a representation which enables them to construct 
solutions in a dynamic programming fashion. We take a general ap- 
proach and relate the construction of such algorithms to the develop- 
ment of algorithms using dynamic programming techniques. Thereby, 
we give general guidelines on how to develop evolutionary algorithms 
that have the additional ability of carrying out dynamic programming 
steps. Finally, we show that for a wide class of the so-called DP- 
benevolent problems (which are known to admit FPTAS) there exists 
a fully polynomial-time randomized approximation scheme based on 
an evolutionary algorithm. 

1 Introduction 

Evolutionary algorithms (EAs) [H] have been shown to be successful for a 
wide range of optimization problems. While these algorithms work well for 
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many optimization problems in practice, a satisfying and rigorous mathe- 
matical understanding of their performance is an important challenge in the 
area of evolutionary computing [1]. 

Interesting results on the runtime behaviour of evolutionary algorithms 
have been obtained for a wide range of combinatorial optimization problems 
(see [31] for a comprehensive presentation). This includes well-known prob- 
lems such as sorting and shortest paths [3S], spanning trees [231, maximum 
matchings [19], and minimum cuts [30| l3T] . There are also some results 
on evolutionary algorithms acting as approximation algorithms for NP-hard 
problems like partition [12], covering [T7], and multi-objective shortest path 
[22l [32] problems. But a general theoretical explanation of the behavior of 
evolutionary algorithms is still missing. The first step in this direction is 
taken in [37], where the authors show for an important subclass of opti- 
mization problems that evolutionary algorithms permit optimal solutions in 
polynomial time. 

1.1 Main Contributions 

The aim of this paper is to make another contribution to the theoretical un- 
derstanding of evolutionary algorithms for combinatorial optimization prob- 
lems. We focus on the question how to represent possible solutions such that 
the search process becomes provably efficient. When designing an evolution- 
ary algorithm for a given problem, a key question is how to choose a good 
representation of possible solutions. This problem has been extensively stud- 
ied in the literature on evolutionary algorithms [38]; for example there are 
different representations for the well-known traveling salesman problem (see 
e.g. Michalewicz [27]) or NP-hard spanning tree problems (see e.g. Raidl 
and Julstrom [3S]). 

Each of these representations induces a different neighborhood of a par- 
ticular solution, and variation operators such as mutation and crossover have 
to be adjusted to the considered representation. Usually, such representa- 
tions either lead directly to feasible solutions for the problem to be optimized 
or the search process is guided towards vahd solutions by using some penalty 
functions. Here, the representation of possible solutions in combination with 
some suitable variation operators may be crucial for the success of the algo- 
rithm. 

Recently, it has been proven for various combinatorial optimization prob- 
lems that they can be solved by evolutionary algorithms in reasonable time 
using a suitable representation together with mutation operators adjusted to 
the given problem. Examples for this approach are the single source shortest 
path problem [39], all-pairs shortest path problem [9], multi-objective short- 
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est path problem [22] , the travelhng salesman problem [H] and the knapsack 
problem [TJ] . The representations used in these papers are different from the 
general encodings working with binary strings as considered earlier in theo- 
retical works on the runtime behavior of evolutionary algorithms. Instead, 
the chosen representations reflect some properties of partial solutions of the 
problem at hand that allow to obtain solutions that can be extended to op- 
timal ones for the considered problem. To obtain such partial solutions the 
algorithms make use of certain diversity mechanisms allowing the algorithms 
to proceed in a dynamic programming way. 

Note that the problem-solving capability of classical genetic algorithms 
is sometimes explained using the building block hypothesis [20], which also 
involves extension of partial solutions to the optimal ones. A relation of the 
mentioned above EAs to dynamic programming, however, allows to obtain 
more specific results in terms of average optimization time. 

Dynamic programming (DP) [3] is a well-known algorithmic technique 
that helps to tackle a wide range of problems. A general framework for 
dynamic programming has been considered by e. g. Woeginger [13] and 
Klotzler [2S]. The technique allows to compute an optimal solution for the 
problem by extending partial solutions to an optimal one. 

An important common feature of the evolutionary algorithms [9l [13], [151 
[22] [391 [H] is that each of them is based on a suitable multi-objective for- 
mulation of the given problem. The schemes of these EAs and solution 
representations are different, however. 

The algorithms proposed in [39] and [15] are generalizations of the well- 
known (l-l-l)-EA (see e.g. [4]) to the multi-objective case and they are based 
on a different representation of solutions than the one used in our paper. 

The (/i -|- 1)-EA in [H] employs a large population of individuals, where 
each individual encodes just one partial solution. In [HI US] it was shown that, 
for the all-pairs shortest path problem on an n-vertex graph, application of 
a suitable crossover operator can provably reduce the optimization time of 
the EA by a factor of almost rr"^^. 

A special case of the DP-based evolutionary algorithm proposed in the 
present paper can be found, e.g. in [22] ■ Both algorithms employ large 
populations of individuals where an individual encodes a partial solution. 
The outline of these algorithms is similar to that of the SEMO algorithm [26]. 

Each gene in our problem representation defines one of the DP transition 
mappings, and a composition of these mappings yields the DP state rep- 
resented by the individual. The proposed EA utilizes a mutation operator 
which is a special case of point mutation, where the gene subject to change is 
not chosen randomly as usual, but selected as the first gene which has never 
been mutated so far (see Section 13.21 for details and links to the biological 
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systems) . 

The goal of the paper is to relate the above mentioned multi-objective 
evolutionary approaches to dynamic programming and give a general setup 
for evolutionary algorithms that are provably able to solve problems having 
a dynamic programming formulation. In particular, we show that in many 
cases a problem that can be solved by dynamic programming in time T has an 
evolutionary algorithm which solves it in expected time 0(T ■ n ■ log {\DP\)) 
with n being the number of phases and \DP\ being the number of states 
produced at the completion of dynamic programming. 

The obtained results are not aimed at the development of faster solution 
methods for the combinatorial optimization problems (to construct an EA 
in our framework, one has to know enough about the problem so that the 
traditional DP algorithm could be applied and this algorithm would be more 
efficient). Instead, we aim at characterizing the area where evolutionary 
algorithms can work efficiently and study the conditions that ensure this. To 
put it informally, our results imply that a class of problems that is easy for 
the DP algorithm is also easy for a suitable EA for most of the reasonable 
meanings of the term "easy" (solvable in polynomial or pseudo-polynomial 
running time or admitting FPTAS). 

1.2 Organization 

The rest of the paper is organized as follows. In Section [21 we introduce 
a general dynamic programming formulation and the kind of problems that 
we want to tackle. This dynamic programming approach is transferred into 
an evolutionary algorithm framework in Section [31 Here we also show how 
to obtain evolutionary algorithms carrying out dynamic programming for 
some well-known combinatorial optimization problems. In Section [H we 
consider a wide class of the DP-benevolent problems which are known to 
have fully polynomial-time approximation schemes based on dynamic pro- 
gramming [l3]. We show that for the problems of this class there exists a 
fully polynomial-time randomized approximation scheme based on an evolu- 
tionary algorithm. Finally, we finish with some conclusions. 

The main results of Sections [H and [3] originally were sketched in our ex- 
tended abstract |8j, while the main result of Section [H was published in Rus- 
sian in [16] . Additionally to refined presentation of results |8l [16] , the present 
paper contains a DP-based EA with a strengthened runtime bound for the 
case of DP algorithm with homogeneous transition functions (applicable e.g. 
to the shortest path problems). 
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2 Dynamic Programming 



Dynamic programming is a general design paradigm for algorithms. The 
basic idea is to divide a problem into subproblems of the same type, and 
to construct a solution for the whole problem using the solutions for the 
subproblems. Dynamic programming has been proven to be effective for 
many single-objective as well as multi-objective optimization problems. It 
is even the most efficient approach known for solution of some problems in 
scheduling [351 US], bioinformatics routing (see e.g. Chapters 24, 25) 
and other areas. 

In this section, we will assume that an original optimization problem 11 
(single-objective or multi-objective) may be transformed into a multi-objective 
optimization problem P of a special type. The general scheme of dynamic 
programming will be presented and studied here in terms of the problem P. 
Several examples of a transformation from 11 to P are provided at the end 
of the section. 

2.1 Mult i- Objective Optimization Problem 

Let us consider a multi-objective optimization problem P which will be well 
suited for application of the DP algorithm in some sense, as shown below. 
Suppose, there are d G N objectives that have to be optimized in P. An 
instance of problem P is defined by a quadruple {d, g,S,V). Here g:S 
(M U {cxo})'^ is called the objective function, S is called the search space, and 
g{S) C (MU {oo})'^ is the objective space. P C 5 is a set of feasible solutions. 

We introduce the following partial order to define the goal in multi- 
objective optimization formally. Throughout this paper, ^ denotes Pareto 
dominance where 

{yi,...,yd) ^ {y'i,...,y'd) 

iff yi > y'i for all i for minimization criteria Qi and yi < y[ for maximization 
criteria gi. In the following, we use the notation y' -< y as an abbreviation 
for y' ^ y and y :^ y'. The Pareto front is the subset of g{'D) that consists 
of all maximal elements of g{T>) with respect to ^. The goal is to deter- 
mine a Pareto- optimal set, that is, a minimal by inclusion subset of feasible 
solutions V that is mapped on the Pareto front. 

2.2 Framework for Dynamic Programs 

Consider a DP algorithm for a problem P, working through a number of 
iterations called phases. In each phase the DP algorithm constructs and 
stores some states belonging to S. By saying that DP algorithm computes 
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a Pareto-optimal set for the problem P we mean that after completion of 
the DP algorithm, the set of all DP states produced at the final phase is a 
Pareto-optimal set for P. 

Application of the DP approach to many multi-objective and single- 
objective optimization problems can be viewed as a transformation of a 
given problem 11 to some problem P: a DP algorithm is applied to com- 
pute a Pareto-optimal set for P and this set is efficiently transformed into a 
solution to the given single- or multi-objective problem. 

In what follows, we consider only those DP algorithms where the states of 
the current phase are computed by means of transition functions, each such 
function depending on the input parameters of problem P and taking as an 
argument some state produced at the previous phase. 

Let us start the formal definition of the DP algorithm from a simplified 
version. Suppose that the simplified DP algorithm works in n phases, such 
that in the i-th phase a set iSj C 5 of states is created. We use n finite 
sets J-'i of state transition functions F:S^S' to describe the DP algorithm. 
Here S' is an extension of space S. A mapping F can produce elements 
F{S) G S'\S that do not belong to a search space. To discard such elements 
at phase i, z = 1, . . . , ra, a consistency function Hi is used. Hi. S' — )■ M, such 
that S* e iS iff Hi{S) < 0. We assume that the number n, the functions Hi 
and the sets of functions J^i depend on the input instance of problem P. 

The simplified DP algorithm proceeds as follows. In the initialization 
phase, the state space Sq is initialized with a finite subset of S. In the i-th 
phase, the state space Si is computed using the state space Si-i according 
to 

Si = {FiS) I 5 G A F G A H,{F{S)) < 0}. (1) 

In the process, the consistency functions Hi serve to keep the infeasible ele- 
ments emerging in phase i from being included into the current state space Si. 
(Note that after completion of phase n of the simplified DP algorithm, the 
set Sn may contain some states whose objective values are Pareto-dominated 
by those of other states from iS„.) 

To delete the states with Pareto-dominated objective values and to im- 
prove the runtime of the simplified DP algorithm defined by ([T]), most of the 
practical DP algorithms utilize the Bellman principle (see e. g. [5]) or its vari- 
ations so as to dismiss unpromising states without affecting the optimality of 
the final set of solutions. A formulation of the Bellman principle in terms of 
recurrence ([1]) for the single-objective problems can be found in|Xl Sufficient 
conditions for application of the Bellman principle in the single-objective case 
were formulated in [2B]- In the multi-objective case the Bellman principle is 
not used, but the unpromising states may be excluded by means of an ap- 
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propriate dominance relation on the set of states. Originally such dominance 
relations were introduced by R. Klotzler [23]. In this paper, we employ a 
similar approach, motivated by |43] . 

Let us consider a partial quasi-order (i. e. a reflexive and transitive rela- 
tion) ^dom deflned on S so that S ^dom S' iff g{S) ^ g{S'). We will say that 
state S is dominated by state S' iff S ^dom S'. If S* G T C 5 is such that 
no S" G T exists satisfying S ^dom S', then S will be called non- dominated 
in T. 

As we will see, under the following two conditions the relation ^dom is 
helpful to dismiss unpromising states in the DP algorithm. 

The flrst condition C.l guarantees that the dominance relation between 
two states transfers from one round to the next: 

Condition C.l. For any 5*, 5" G = l,...,n, if S ^dom S' then 

FiS) ^dom F{S') for all F G Ti. 

The second condition C.2 expresses that infeasible states cannot dominate 
feasible states: 

Condition C.2. For any S,S' G S', if S ^dom S' and Hi{S) < then 
HiiS') < 0. 

Consider a subset Si of S. We call % Si a dominating subset of Si 
with respect to ^dom iff for any state S E Si there is a state S" G 71 with 
•S* ^dom S'. Let us use the notation M {Si, ^dom) to denote the set of all 
dominating subsets of Si which are minimal by inclusion. 

The following proposition indicates that under conditions C.l and C.2 it 
is sufficient to keep a dominating subset of states constructed in each phase i, 
rather than the full subset Si. 

Proposition 1. Suppose the simplified DP algorithm is defined by condi- 
tions C.l and C.2 hold and the dominating sets 71, i = 1, . . . ,n are computed 
so that % G M{So, ^dom), 

71 G M{{F{S) \SeTi-iAFeTiA H,{F{S)) < 0}, ^dom). (2) 

Then for any state S* E Si, i = 0, . . . ,n, there exists S E % such that 

S* ^dom S . 

Proof. The proof is by induction on i. For i = the statement holds by 
assumption % G M{So, ^dom)- 

By ([I]), a state S* G Si can be expressed as S* = F*{S'), so that 
Hi{S*) <0,F* G J-'i and S' G But by induction hypothesis, there 

exists a state G 7i_i such that S' ^dom S^- Now conditions C.l and C.2 
imply that S* = F*{S') ^dom F*{S^) and F*{S^) G {F{S) | 5 G 71-1 A F G 
J^iAHi{F{S)) < 0}. Hence, by ([2]) we conclude that there exists S E Ti such 



7 



Algorithm 1 Dynamic Program for P 

1: To^fH 

2: for S eSodo 

3: if $S' eTq-.S ^dom S' then 

4: To ^ (To \ {S' eTolS' ^dom S}) u {5} 

5: end if 

6: end for 

7: for 2 = 1 to n do 

8: 71^0 

9: for S G 71-1 and F G -Fj do 

10: if Hi{F{S)) < and ^5' G 71: F{S) ^dom 5' then 

11: 7- ^ (7; \ G 71 I 5' ^dom F{S)}) U 

12: end if 

13: end for 

14: end for 

15: return Tn 



that S* ^dom ^dom 5. □ 

In view of definition of ^dom, if the conditions of Proposition[T]are satisfied 
and the Pareto front of g is contained in g{Sn), then this Pareto front is also 
contained in g{Tn)- 

Proposition 2. If the conditions of Proposition^ are satisfied, then the size 
of each Ti, i = 0, . . . ,n, is uniquely determined. 

Indeed, consider the set of maximal elements of Si with respect to ^dom- 
Define the equivalence classes of this set with respect to the equivalence 
relation x = y iS x ^dom y and y ^dom The size of a minimal subset M of 
the set of maximal elements, which dominates all elements of Si, is unique 
since such M contains one representative element from each equivalence class. 

□ 

A computation satisfying ([2]) can be expressed in an algorithmic form as 
presented in Algorithm [TJ It is easy to see that when a subset % is completed 
in Lines 8-13, condition holds. 

The runtime of a DP algorithm depends on the computation times for 
the state transition functions F G J-^, for the consistency functions Hi, for 
checking the dominance and manipulations with the sets of states. Let ^i? be 
an upper bound on computation time for a transition function F and let 6-^ 
be an upper bound for computation time of any function Hi, i = 1, . . . ,n. 
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Sometimes it will be appropriate to use the average computation time for 
the state transition functions at phase i, i = 1, . . . ,n: 9jr- = J^PeTi ^F/\J^i\- 
In Algorithm [H verification of condition 

$S'e%:F{S) ^dom^' (3) 

in Line 10 and execution of Line 11 may be implemented using similar 
problem-specific data structures. To take this into account, we will denote 
by 6^ an upper bound applicable both for the time to verify ([3]) and for the 
time to execute Line 11. 

The body (lines 10-12) of the main loop (Lines 7-14) in Algorithm [T] is 
executed Yl^=i I'^i-il ' l-^il times. 

To simplify the subsequent analysis let us assume that in the if-statement 
at line 10, the condition ([3j) is always checked. We denote the computation 
time for initializing % with 6*1111 (Lines 1-6) and the computation time for 
presenting the result with 6'out (Line 15), which leads to an overall runtime 

O ^^ini + J2 ■ I'^-ll ■ (^-^^ + ^« + + j • (4) 

In many applications of the DP, the computation time for the state tran- 
sition functions and the consistency functions are constant. Besides that, 
the partial quasi-order ^dom is often just a product of linear orders and it is 
sufficient to allocate one element in a memory array to store one (best found) 
element for each of the linear orders. This data structure usually allows to 
verify ([3]) and to execute Line 11 in constant time (see the examples in Sub- 
section [HIS])- In the cases mentioned above, the values 9p, Oy^ and 9^ can 
be chosen equal to the corresponding computation times and the overall DP 
algorithm runtime in (jlj) can be expressed with symbol G(-) instead of O(-). 

Note that the runtime of the DP algorithm is polynomially bounded in 
the input length of problem P if 6'ini, n, 6jr-, 6'-^, 9^, 6'out, as well as |7I| and 
for i = 0, ... ,n — 1, are polynomially bounded in the input length. 
Here and below, we say that a value (e.g. the running time) is polynomially 
bounded in the input length, meaning that there exists a polynomial function 
of the input length, which bounds the value from above. 

2.3 Applications of the general DP scheme 

In this subsection, we point out how the general DP framework presented 
above is applied to some classical combinatorial optimization problems. The 
approach followed here is to describe the appropriate problem P and the 



9 



components of a dynamic programming algorithm for the solution of a spe- 
cific problem 11. Most of the following examples have been inspired by the 
previous works [HI |39l E] . Note that ^dom will be a product of linear orders 
in each of these examples. In what follows id denotes the identical mapping. 

Traveling Salesman Problem Let us first consider the traveling sales- 
man problem (TSP) as a prominent NP-hard example. The input for the 
TSP consists of a complete graph Q = (V, E) with a set of nodes V = 
{1,2, ... ,n} and non-negative edge weights w: E ^ Mg. It is required to 
find a permutation of all nodes (fi, . . . ,Vn), such that the TSP tour length 
X]r=2 ^(^i-i' + '^iyn.vij is minimized. Without loss of generality we can 
assume that fi = 1, that is, the TSP tour starts in the fixed vertex 1. 

The search space S for problem P corresponding to the dynamic program- 
ming algorithm of Held and Karp [21] consists of all paths S = [vi, . . . ,Vi),Vi = 
1 of i = 1, . . . ,n nodes. S' is the extended search space of all sequences of 
nodes up to length n (the same node may occur more than once). Given 
M C V\{1} and k e M, let Ti{k, M) denote the set of all paths of |M| + 1 
vertices starting in vertex 1 then running over all nodes from M and end- 
ing in vertex k. Let the vector objective function (7 : 5 — > (M U {00})'^, 
d = {n — 1)2^""^^ have components gkAiiS) for all M C V^\{1}, k G M, equal 
to the length of path 5* iff S* G TT{k,M). For all other S ^ n{k,M) assume 
gkAiiS) = 00. The set of feasible solutions is V = U^^2^(^5 since in 

the TSP we seek a tour of length n. 

Sq consists of a single element Vi. The set J-", for all i consists of n — 1 
functions F^-.S S' that add vertex v G to the end of the given path. 

For invalid states S & S', which are characterized by not being Hamiltonian 
paths on their vertex sets, the mapping Hi{S) computes 1 and otherwise. 

In view of the definition of objective g, the dominance relation is formu- 
lated as follows. S ^dom S' if and only if S and S' are Hamiltonian paths on 
the same ground set with the same end vertex k and path S' is not longer 
than S. States from different sets 7i{k,M) are not comparable. Conditions 
C.l and C.2 are verified straightforwardly. 

Substituting these components into Algorithm [H we get almost the whole 
well-known dynamic programming algorithm of Held and Karp [21] , except 
for the last step where the optimal tour is constructed from the optimal 
Hamiltonian paths. 

Algorithm [1] initializes the states of the dynamic program with paths 
(1, f ) for all G \ {1}. In each subsequent iteration i, the algorithm takes 
each partial solution S obtained in the preceding iteration and checks for 
every application of the state transition function F{S) with F & J^i whether 
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Hi{F{S)) is a feasible partial solution that is non-dominated in %. If so, then 
F{S) is added to the set % of new partial solutions by replacing dominated 
partial solutions S' defined on the same ground set with the same end vertex 
of the Hamiltonian path. 

What remains to do after completion of the DP algorithm with Pareto- 
optimal set is to output the Pareto-optimal solution minimizing the criterion 
9k,v\{i}{S) + w{k,l),k G which is now easy to find. Here using 

appropriate data structures one gets 9f = ©(1), = ©(l); = ©(1) 
and \Ti\ = \J-'i\ = n — 1 for all i = l,...,n, thus the observation 

following (jlj) leads to the time complexity bound ©(^^2"^). 

Knapsack Problem Another well-known NP-hard combinatorial optimiza- 
tion problem that can be solved by dynamic programming is the knapsack 
problem. The input for the knapsack problem consists of n items where each 
item i has an associated integer weight Wi > and profit pi > 0, 1 < i < n. 
Additionally a weight bound W is given. The goal is to determine an item 
selection K C {1, . . . ,n} that maximizes the profit Ylii<^KPi^ subject to the 
condition Xligx"^* — 

We fit the problem into the above framework assuming that each state 
S = (si,S2) G Si, i = l,...,n, encodes a partial solution for the first i 
items, where coordinate si stands for the weight of a partial solution and S2 
is its profit. The initial set So consists of a single element (0, 0) encoding a 
selection of no items. 

The pseudo-Boolean vector function g:S ^ M.^ defines W criteria 

oieJise • " = «).-.^''. (5) 

that have to be maximized. This implies the dominance relation :<dom such 
that 5* ^dom S' iff si = s[ and S2 < where 5* = (si, S2), S' = {s[, S2). 

The set J-i consists of two functions: id and Fi{si, S2) = (si +Wi, S2 +Pi)- 
Here Fi corresponds to adding the z-th item to the partial solution, and id 
corresponds to skipping this item. A new state S = (si, S2) is accepted if it 
does not violate the weight limit, i. e. Hi{S) < 0, where Hi{S) = Si — W. 

The conditions C.l and C.2 are straightforwardly verified. To obtain an 
optimal solution for the knapsack problem it suffices to select the Pareto- 
optimal state with a maximal component S2 from Sn- 

To reduce the comparison time 6^ we can store the states of the DP in 
a {W X n)-matrix. An element in row w, w = 1, . . . , W, and column i, i = 
1, . . . , n, holds the best value S2 obtained so far on states S = (si, S2) G % 
with si = w. Then 6^ is a constant and the worst-case runtime of the 
explained DP algorithm is 0{n ■ W) since Y17=i I'^i-^l — 
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Single Source Shortest Path Problem A classical problem that also fits 
into the DP framework is the single source shortest path problem (SSSP). 
Given an undirected connected graph Q = (V, E), \V\ = n and positive edge 
weights w: E ^ M"*", the task is to find shortest paths from a selected source 
vertex s G to all other vertices. 

The search space iS is a set of all paths in Q with an end-point s. The set 
of feasible solutions V is equal to S. 

Since adding a vertex to a path may result in a sequence of vertices that 
do not constitute a path in Q, we extend the search space to the set S' 
of all sequences of vertices of length at most n with an end-point s. The 
set Sq of initial solutions is just a single vertex s. Now for all i, we define 
J^i := {Fy \ V E V} U {id}, where F^:S — t- S' is the mapping adding the 
vertex f to a sequence of vertices. Hi{S) = —1 if is a path in Q with an 
end-point s, and 1 if not. 

Let the vector objective function g have d = n components gv{S) for all 
V E V, equal to the length of path S iff S connects s to v, otherwise assume 
gv{S) = CO. This implies that S ^dom S' if and only if the paths S and S' 
connect s to the same vertex and S' is not longer than S. 

The resulting DP algorithm has Op = 0(1), 0^ = 6(1), = 9(1) and 
\Ti-i\ = 6(n), \J^i\ = 9(n) for all i = l,...,n, thus (j4]) gives the time 
complexity bound 0{n^). The well-known Dijkstra's algorithm has 0(?t,^) 
time bound, but in that algorithm only one transition mapping is applied in 
each phase (attaching the closest vertex to the set of already reached ones), 
and such a problem-specific DP scheme is not considered here. 

All-Pairs Shortest Path Problem Finally, let us consider the all-pairs 
shortest path (APSP) problem, which has the same input as the SSSP, except 
that no source vertex is given, and the goal is to find for each pair {u, v) of 
vertices a shortest path connecting them. 

A basic observation is that sub-paths of shortest paths are shortest paths 
again. Hence a shortest path connecting u and v can be obtained from 
appending the edge (x, f), where a; is a neighbor of f, to a shortest path 
from u to X. This allows a very natural DP formulation as described for 
problem P. 

For the APSP, the search space S naturally is the set of all paths in ^, 
and the set V of feasible solutions consists of collections of paths, where for 
each pair of vertices there is one path connecting them. 

We model paths via finite sequences of vertices, and do not allow cycles. 
Since adding a vertex to a path may create a sequence of vertices which does 
not correspond to a path in let us extend this search space to the set S' of 
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all sequences of vertices of length at most n. The set Sq of initial solutions is 
the set of all paths of length 0, that is, of all sequences consisting of a single 
vertex. Now for all i, we define J-i := {F^ \ v G V} U {id}, where F^-.S' S' 
is the mapping adding the vertex f to a sequence of vertices. To exclude 
invalid solutions, let us define Hi{S) to be —1 if S" is a path in Q, and 1 if 
not. 

It remains to define when one state dominates another. Let vTjj denote 
the set of all paths starting in vertex i and ending in vertex j. Let the vector 
objective function g : S ^ {RU {oo})'^, d = have components gij{S) for 
all i,j e V, equal to the length of path S iS S E T^ij. For all other S ^ nij 
assume gij{S) = oo. This implies that 5* ^dom S' if and only if the paths S 
and S' connect the same two vertices and 5" is not longer than S. 

Since the length of the path arising from extending an existing path by 
an edge depends monotonically on the length of the existing path, conditions 
C.l and C.2 hold. So, in view of Proposition [H any set 71 contains a path 
for each pair of vertices (and only one such path). Thus, 7^ is a subset of V 
and contains a shortest path for any pair of vertices. 

The resulting algorithm following the dynamic programming approach 
now does the following. It starts with all paths of length zero as solution 
set Sq. It then repeats n times the following. For each path in the solution 
set and each vertex, it appends the vertex to the path. If the resulting path 
dominates an existing solution with the same end vertices, it replaces the 
latter. Here Op = 9(1), On = 0(1), 6-, = 6(1) and |7I| = ©(n^), | = 0{n) 
for alH = 1, . . . ,n, thus (j4j) gives the time complexity bound O(n^). Note 
that the well-known Floyd- Warshall algorithm (see e.g. |7j. Chapter 25) has 
O(n^) time bound, but in that algorithm each transition mapping combines 
two states (paths), and such an option is not considered in this paper. 

3 Evolutionary Algorithms 

In the following, we show how results of dynamic programming can be at- 
tained by evolutionary algorithms. To this aim, we state a general formula- 
tion of such an evolutionary algorithm and then describe how the different 
components have to be designed. 

3.1 Framework for Evolutionary Algorithms 

An evolutionary algorithm consists of different generic modules, which have 
to be made precise by the user to best fit to the problem. Experimental prac- 
tice, but also some theoretical work (see e.g. [IHl [HI [121 129]), demonstrate 
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Algorithm 2 Evolutionary Algorithm for P 



for I eVo do 

if $r eV:I ^EA /' then 

v^{v\{r ev\r /}) u {/} 

end if 
end for 
loop 

/ <- mut(sel(P)) 

if $r eV:I ^EA /' then 

v^iv\{r ev\r ^EA /}) u {/} 

end if 
end loop 

return {outEA(/) \ I = (i, S) E V, S E V} 



that the right choice of representation, variation operators, and selection 
method is crucial for the success of such algorithms. 

We assume again that an instance of problem P is given by a multi- 
objective function g that has to be optimized. We consider simple evolution- 
ary algorithms that consist of the following components. 

We use iS^A := {0, . . . , n} x S' as the phenotype space and call its elements 
individuals. The algorithm (see Algorithmic]) starts with an initial population 
of individuals Vo- During the optimization the evolutionary algorithm uses 
a selection operator sel(-) and a mutation operator mut(-) to create new 
individuals. The dimensional objective function together with a partial 
order ^ on M'^ induce a partial quasi-order ^ea on the phenotype space, 
which guides the search. After the termination of the EA, an output function 
outEA(') is utilized to map the individuals in the last population to search 
points from the DP search space. 

3.2 Defining the Modules 

We now consider how the different modules of the evolutionary algorithm 
have to be implemented so that it can carry out dynamic programming. To 
do this, we relate the modules to the different components of a DP algorithm. 
Consider a problem P given by a set of feasible solutions V and a multi- 
objective function g that can be solved by a dynamic programming approach. 
The EA works with the following setting. 

The initial population is Vq = {0} x Sq where So is the initial state space 
of the DP algorithm. The selection operator sel(-) chooses an individual 
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I the following way. First it chooses i < n — 1 uniformly from the set 
of phases which are represented in the current population i. e. from the set 
{k : k < n — 1, 3(A;, S) G V}. After this, selection chooses / uniformly among 
the individuals of the form (i, S) in the current population. 

For an individual {i,S), the mutation operator mut(-) chooses a state 
transition function F G J-i+i uniformly at random and sets mut((i, S*)) = 
{i + l,F{S)). 

We incorporate a partial order ^ea into the EA to guide the search. This 
relation is defined as follows: 

S) ^EA {i\ S') ^{i = i' and S ^dom S') or H,{S) > 0. (6) 

Finally, we utilize the output function outEA((^, 5)) = 5* to remove the ad- 
ditional information at the end of a run of the EA. That is, we remove the 
information that was used to store the number of a certain round of the un- 
derlying dynamic program and transform an individual into a search point 
for the problem P. 

Note that the description of the Algorithm [2] does not employ the notion of 
the fitness function, although an appropriate multi-objective fitness function 
may be defined for compatibility with the standard EA terminology. 

Finally, note that we do not discuss the solutions encoding in our EA 
because it is not essential for the analysis. However, it may be worth men- 
tioning, when the biological analogy is considered. Here each of the genes 
Ai, i = 1, . . . ,n would define the DP transition mapping from a set J-^, and 
a composition of these mappings would yield the DP state represented by 
the individual. One of the possible options of each gene is "undefined" , and 
the mutation operator modifies the first gene which is still "undefined" in 
the parent individual. A discussion of genetic mechanisms corresponding to 
the proposed mutation in a biological system is provided in |Bl 

3.3 Runtime of the Evolutionary Algorithm 

Our goal is to show that the evolutionary algorithm solves the problem P 
efficiently if the dynamic programming approach does. To measure the time 
the evolutionary algorithm needs to compute a Pareto-optimal set for prob- 
lem P, one would analyze the expected number of fitness evaluations to come 
up with a Pareto-optimal set, when it is non-empty. This is also called the 
expected optimization time, which is a common measure for analyzing the 
runtime behavior of evolutionary algorithms. The proposed EA does not use 
the multi-objective fitness function explicitly, but given enough memory, it 
may be implemented so that every individual constructed and evaluated in 
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Lines 3 and 4 or in Lines 8-11 requires at most one evaluation of the objective 
function g. Thus, we can define the optimization time for Algorithm [2] as 
|iSo| plus the number of iterations of the main loop (Lines 8-11) required to 
come up with a Pareto-optimal set. Analogous parameter of a DP algorithm 
is the number of states computed during its execution. 

The next theorem relates the expected optimization time of the EA to 
the number of states computed during the execution of the corresponding 
DP algorithm. In what follows it will be convenient to denote the cardinality 
of the set of states produced after completion of the DP algorithm by \DP\, 
i.e. \DP\ := Ym=o I'^I- Note that \DP\ is a well-defined value since the sizes 
\Ti\ are unique according to Proposition [2j 

Theorem 1. Let a DP he defined as in AlgorithmU\ and an EA defined as 
in Algorithmic with :<ea relation defined by Then the number of states 
computed during the execution of the DP algorithm is \So\ +Yl^=i ' I'^i-il; 
and the EA has an expected optimization time of 

O (^So\ + n ■ log \DP\ ■ I'^l ■ l-^^+il j • 

Proof. Estimation of the number of states computed during the execu- 
tion of the DP algorithm is straightforward. 

Assume that the optimization process works in stages 1 < i < n, whereas 
stage i + 1 starts after the stage i has been finished. We define that a stage i 
finishes when for every state S E % there exists an individual {i, S') G V 
with S' dominating S. Here and below in this proof, by 71, i = 0, . . . ,n, 
we denote the corresponding sets computed in Algorithm [H Note that after 
completion of a stage i, the subset of individuals of a form (i, S) in population 
V does not change in the subsequent iterations of the EA. Let 1~' denote the 
set of states of these individuals after completion of stage i, i = 0, . . . ,n. By 
the definition of Algorithm [21 the sequence i = 0, . . . ,n satisfies (j2]), and 
therefore = |7i|, i = 0, . . . ,n in view of Proposition [2l 

Let be the random variable denoting the number of iterations since 
stage z — 1 is finished, until stage i is completed. Then the expected opti- 
mization time is given by |iSo| + E[^] with ^ = + . . . + ^n- 

Any state S G 7i+i is computed in Algorithm [1] by means of some function 
F G J^i+i, when it is applied to some state S E %■ Thus, in stage i -|- 1 of 
the EA during mutation the same transition function F may be applied to 
some individual /' = {i,S'), such that S ^dom S'. After this mutation, in 
view of conditions C.l and C.2, the population V will contain an individual 

J" = (Z + 1, S") with S" such that S ^dom F{S') ^dom S". 
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Consider any iteration of the EA at stage i + Let t denote the number of 
such states from TI+i that are already dominated by a state of some individual 
in V. Then there should be |7i+i| — t new individuals of the form {i + 1, S) 
to be added into V to complete stage i + 1 (recall that |T'i+i| = |7i+i|). The 
probability to produce an individual {i, S') where S' dominates a previously 
non-dominated state from 7i+i is no less than (|7i+i| — t)/{n\Ti\ ■ 
with an expected waiting time of at most (n|7i| • |~^j+i|)/(|7I+i| — t) for this 
geometrically distributed variable. The expected waiting time to finish stage 
z + 1 is thus bounded by 

E[^i+i] < ^ = n\%\ ■ \J^i+i\ ■ H|7-+i|, 

t=i 

with 'Hk being the fc-th harmonic number, Tik := Ylt=i 1- 
This leads to an overall expected number of iterations 

n— 1 n— 1 

E[^] < 5^n|7;| ■ I J-i+il ■ H|7-+,| < n{\n \DP\ + 1) . ^ IT"] • | 

1=0 i=0 

□ 

A similar inspection as in S ubsect ion 12 . 2 1 re veals that the expected runtime 
of the EA is 

n-l 

0(^ini + nlog \DP\ ■ J2 (l-^^+il ■ I'TII ■ (^^: + 0h + ^d)) + ^out) , 

assuming the individuals of the population are stored in + 1 disjoint sets 
according to the first coordinate i. 

As noted in Subsection 12.21 if the computation times for functions F, 
Hi and dominance checking as well as execution time for Line 11 in 
Algorithm [T] are constant, then 6p, 6'-^ and 6^ can be chosen equal to the 
corresponding computation times. In such problem that is solved 

by dynamic programming Algorithm [1] in time T, will be solved by the EA 
defined as in Algorithm [2] in expected time 0(Tn log \ DP\). 

3.4 Homogeneous transitions 

Some DP algorithms, like the ones for the APSP and SSSP problems, have 
a specific structure which may be exploited in the EA. In this subsection we 
consider the case of homogeneous transition functions where T\ = ...=Tn 
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and Hi = . . . = Hn- To simplify the notation in this case we will assume 
T\ = T and Hi{S) = H{S). Additionally, we suppose that the identical 
mapping belongs to J-". 

The formulated assumptions imply that once some state S is obtained 
in the DP algorithm, it will be copied from one phase to another, unless 
some other state will dominate it. Note also that it does not matter at 
what particular phase a state has been obtained - the transition functions 
will produce the same images of this state. These observations motivate 
a modification of the partial order ^ea, neglecting the phase number in 
comparison of individuals: 

{i, S) ^EA it', S') ^ S ^dom S' or H^iS) > 0. (7) 

In fact, now we can skip the index i in individuals [i, S) of the EA, so in 
this subsection the terms "state" and "individual" are synonyms and the 
phase number i is suppressed in the notation of individuals. As the following 
theorem shows, wider sets of comparable individuals in this special case allow 
to reduce the population size and thus improve the performance of the EA. 
Let us consider the width VFdom of partial order ^dom, i- e. the maximum size 
of a set of pairwise incomparable elements. 

Theorem 2. // the transition functions are homogeneous and id G then 
the EA defined as in Algorithmic with the modified ^ea relation Ij^ has an 
expected optimization time of 0{\So\ + Wdom log(Wdom) ■ '"■l-^D- 

Proof. The analysis is similar to the proof of Theorem [TJ Note that now 
the size of population V does not exceed Wdom- We assume that 1^1 = Wdom 
right from the start. 

Let 71 be the same as in phase i of the DP algorithm, i = 0, . . . ,n. Sup- 
pose again that the optimization process works in stages 1 < i < n, whereas 
stage i is assumed to be finished when for every S E the population V 
contains an individual 5" such that S ^dom S' . 

Let be the number of iterations since stage i — 1 is finished, until stage 
i is completed. Then the expected optimization time is given by \Sq\ + E[^] 

with e = 6 + • • • + ^n- 

Any state 5* G 7i+i is computed in the DP algorithm by means of some 
function F G J^j+i, when it is applied to some state S" G 71. Thus, in stage i+1 
of the EA during mutation the same transition function F may be applied 
to some individual I' = S', such that S ^dom S'. After this mutation, in 
view of conditions C.l and C.2, the population V will contain an individual 

/" = S" such that S ^dom F{S') ^dom 

The probability of such a mutation for a particular S G 7i+i is at least 
1/(1 ■ |P|) < 1/(1 -FFdom)- Let t denote the number of states 5 G 7i+i 
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that are already dominated at stage i + 1. Then there are at least |7i+i| — t 
possibilities to add a new individual, which dominates a previously non- 
dominated state from 7i+i. The probability for such a mutation is not less 
than (|7i+i| — t)/(Wdora ■ with an expected waiting time of at most 

(W^dom- \J^i+i\)/ —t) for this geometrically distributed variable. The ex- 
pected waiting time to finish stage i + 1 is thus E'f^j+i] < W^om- '"^iT^+il- 
But |7i+i| < Wdom because the states of TI+i are pairwise incomparable ac- 
cording to Algorithm [H This leads to an overall expected number of itera- 
tions E[^] < W^dom ■ (In(iydom) + 1) ■ Er=o' l-^i+il- n 

3.5 Examples 

Now, we point out how the framework presented in this section can be used 
to construct evolutionary algorithms using the examples from Section O 

Traveling Salesman Problem Due to Theorem [1] the expected optimiza- 
tion time of the evolutionary algorithm based on the DP algorithm of Held 
and Karp presented in Section [2.31 is 0(n^2"). This bound can be further 
improved to 0{n^ ■ 2") for the EA proposed in \41\ . 

Knapsack Problem Consider the DP algorithm presented in Section 12. 3[ 
The expected optimization time of the corresponding EA for the knapsack 
problem is 0(n^ ■ W ■ log {n ■ W)) due to Theorem [TJ 

Single Source Shortest Path Problem Application of Theorem [1] to 
the DP algorithm for SSSP problem from Section [2^31 gives an expected op- 
timization time of O (n^ log (n) ) for Algorithm [2l 

The DP algorithm for SSSP problem has homogeneous transition func- 
tions with with W^om = Thus, the modified EA considered in Theorem [2] 
has the expected optimization time 0{n^logn). This bound can be further 
improved to 0{n^) for the (1+1)-EA |39j . 

All-Pairs Shortest Path Problem Plugging the ideas of the DP algo- 
rithm for APSP problem presented in Section [231 into the framework of Algo- 
rithm[2], we obtain an EA with an expected optimization time of 0{n^ log (n)) 
due to Theorem [H 

It has been noted, however, that the DP algorithm for APSP has homoge- 
neous transition functions, each set J-i contains the identical mapping. Here 
Wdom = n^, thus Theorem [2] implies that the modified EA has the expected 
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optimization time 0{n^ logn). This algorithm can be further improved to an 
EA with optimization time Q{n^) as has been shown in [9]. 

4 Approximation Schemes 

In this section, we demonstrate that for many single-objective discrete op- 
timization problems 11 the above framework can be used to find feasible 
solutions with any desired precision. The supplementary multi-objective 
problem P will be formally introduced for compatibility with the previous 
sections, but it will not play a significant role here. 

Throughout this section we assume that 11 is an NP-optimization prob- 
lem [2], X denotes the input data of an instance of 11, Sol^ is the set of feasible 
solutions, mx : Sol^ — j- Nq is the objective function (here and below No de- 
notes the set of non-negative integers). The optimal value of the objective 
function is OPT{x) = max^^^^^ ^x{y) if H is a maximization problem, or 
OPT{x.) = min^g^^^ fn^iv) in the case of minimization. To simplify presen- 
tation in what follows we assume that Sol^ 7^ 0. 

To formulate the main result of this section let us start with two standard 
definitions [18]. 

A p- approximation algorithm for 11 is an algorithm that for any instance x 
returns a feasible solution whose objective value at most p times deviates 
from OPT{x) (if the instance x is solvable). Such a solution is called p- 
approximate. A fully polynomial time approximation scheme (FPTAS) for a 
problem 11 is a family of (1 + e)- approximation algorithms over all factors 
e > with polynomially bounded running time in problem input size |x| and 
in 1/e. 

In [l3] G. Woeginger proposed a very general FPTAS with an outline 
similar to the DP Algorithm [H except that the comparison of newly gener- 
ated states to the former ones is modified so that the "close" states are not 
kept. This modified algorithm is denoted by DP a in what follows (a detailed 
description of DP a will be given in Subsection 14.11) . 

The state space S and its subsets 71 computed in the DP Algorithm [T] may 
be exponential in problem input size, thus leading to an exponential running 
time of the DP algorithm (this holds e.g. for the Knapsack problem). The 
algorithm DP a, however, iteratively thins out the state space of the dynamic 
program and substitutes the states that are "close" to each other by a single 
representative, thus bringing the size of the subsets 71 down to polynomial. 
This transformation is known as trimming the state space approach. 

In [l3], a list of conditions is presented, that guarantee the existence of an 
FPTAS when there is an exact DP algorithm for a problem. If a problem 11 
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satisfies these conditions, it is called DP-benevolent. This class, in particu- 
lar, contains the knapsack problem and different scheduling problems, e.g. 
minimizing the total weighted job completion time on a constant number 
of parallel machines, minimizing weighted earliness-tardiness about a com- 
mon non- restrictive due date on a single machine, minimizing the weighted 
number of tardy jobs etc. The definition of DP-benevolence is as follows. 

The input data of 11 has to be structured so that x consists of n vectors 
Xi, . . . , Xn G Nq and the components Xu, . . . , Xai of each vector Xi are given 
in binary coding. The dimension a may depend on the specific problem 
input. 

Suppose that for a problem 11 there exists a corresponding multi-objective 
problem P and an exact simplified DP algorithm defined by expression ([T]). 
This algorithm works in n phases and for each i = 1, . . . ,n the set of func- 
tions J^i and the function Hi do not depend on any input vectors other 
than Xi. Besides that, S G S' = Nq, where dimension /3 is fixed for 11 and 
does not depend on a particular input x. The assumption that elements of S' 
are integer vectors will be essential in this section because each component 
of a state will actually be a quantitative parameter and will be subject to 
scaling. It is sometimes possible, however, to move from integer components 
to reals using the approach from [5]. 

The reduction from 11 to P, according to Section |2l implies that the 
Pareto-optimal set of P can be efficiently transformed into a solution to 
the problem 11. Now let us suppose additionally that any 5* G 5„ can be 
mapped to some y{S) G Sol^ and there is a function G : Nq — i- No such that 
m^iy{S)) = G{S). 

The assumption that the simplified DP algorithm described in Section |2] 
provides an exact solution to 11 may be expressed formally: 

(9Pr(x) = mm{G{S) : S G (8) 

if n is a minimization problem, or alternatively 

OPr(x) = max{G{S) : S G S^}, (9) 

if n is a maximization problem. 

The function y{S) is usually computed by means of a standard back- 
tracking procedure (see e.g. [7], Chapter 15). A general description of such 
a procedure is beyond the scope of the paper since the details of reduction 
from problem 11 to P are not considered here. 

Suppose a degree vector D = {di, . . . , dp) G Nq is defined for 11. Then, 
given a real value A > 1 we say that S = (si,...,^/?) is {D, A)-close to 
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S' = {s[,...,s'p),if 

A-'^^S£ < 4 < A'^^s^, i=l,...,p. 

Let us denote by Cq the set of indices I < £ < /3 such that = 0, and let 

£i = {l,...,/3}\£o. 

The main tool to exclude unpromising states in a DP-based FPTAS |13] 
is the quasi-linear order ^qua, which is an extension of a partial order ^dom, 
i. e. if S ^dom 5" then S ^qua S' for any S, S' eNq. For the sake of compat- 
ibility with |l3], we will limit the consideration to the case where ^dom is a 
partial order, rather than a more general partial quasi-order as in Sections |2] 
and [3l This restriction is not significant w. r. t. applications of the frame- 
work, although most likely the results of [13], as well as our results below, 
hold for the partial quasi-orders as well. 

At each phase z, i = 1, . . . , n, in DP a only those states S may be excluded 
that are dominated in terms of ^qua by one of the other obtained states S', 
provided that S' is {D, A)-close to S. 

Note that for any instance x the partial order ^dom on the final sets 
Si, . . . ,Sn niay be represented by a finite number of criteria gi, . . . ^g^ of 
a corresponding instance of the problem P so that the Pareto-dominance 
relation is equivalent to ^dom on this set. 

A problem 11 is called DP-henevolentiihesidi.es C.l and C.2, the following 
conditions C.1',C.2',C.3 and C.4 hold: 

Condition C.l'. For any A > 1, 5, 5' G and F G Ji, i = 1, . . . , n, if S 
is (D, A)-close to S' and S ^q^a S\ then either F{S) ^q^a F{S') and F{S) 
is (D, A)-close to or F{S) ^dom F{S'). 

Condition C.2'. For any A > 1, 5, 5' G and z = if 5 is 

(D, A)-close to S' and S ^q^a S' , then Hi{S') < Hi{S). 

Condition C.3. A value 7 G No exists, depending only on G and D, such 
that for any A > 1 and S, S' eN^, 

(i) if S is {D, A)-close to S' and S ^q^a S', then G{S') < A^G{S) in 
the case of minimization, and A^"'G{S) < G{S') in the case of maximization 
problem, 

(ii) if S ^dom 5", then G{S') < G{S) in the case of minimization, and 
G{S') > G{S) in the case of maximization problem. 
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Condition C.4. 

(i) The functions F G J-i, iJj, z = 1, . . . , and G, as well as the relation 
^qua are computable in time polynomially bounded in the input length. 

(ii) I i = 1, . . . , n is polynomially bounded in input length. 

(iii) (So is computable in time polynomially bounded in input length. 

(iv) A polynomial 7ri(n, logg |x|) exists, such that all coordinates of any 
element S" G iSj, i = 1, . . . , n are integer numbers bounded by e'^^*^"'''"^^ l^l)^ 
Besides that, for all £ G the cardinality of the set of values that such a 
coordinate can take : (si, . . . , s^, . . . , Sfi) G Si}\ is bounded by a polyno- 
mial 7r2(n, log2 |x|). 

Example: knapsack problem We can verify the DP-benevolence con- 
ditions for the knapsack problem as a simple illustrating example. Let the 
problem input, the DP states and the sets of mappings J-^, i = 1, . . . as 
well as functions Hi be defined as in Section 12.31 Besides that, G{S) = S2 
for all S = (si, S2) G Sn and the degree vector is D = (1, 1). 

A proper linear quasi-order :<qua that suits the partial order ^dom defined 
in Section 12.31 for the knapsack problem is not known to us. Instead, we can 
consider the following relations ^qua and ^dom: let S ^qua S' iff Si > s[, 
where S = (si,S2), S' = {s[,S2) and let ^dom be the trivial partial order, 
i. e. S ^dom 5" iff 5 = S'. (For an example of a DP-benevolent problem with 
non-trivial ^dom see the problem of minimizing total late work on a single 
machine [43].) 

The statements in Conditions C.l, C.2, and C.3(ii) are fulfilled since 
^dom is trivial. The function G{si, S2) = S2 satisfies Condition C.3(i), which 
can be verified straightforwardly, assuming 7 = 1. To see that Condition 
C.4 holds, consider a polynomial 7ri(n, log2 |x|) = ln(2l^l), which ensures that 
max{s^ G 5i|£= 1,2, i = l,...,n} < max{^"^^pi, M^} < 2^ = e^iW°g2|xl), 

Conditions C.l' and C.2' hold because the functions Fi, id and Hi at any 
phase i just sum the arguments with given non-negative constants. Indeed, 
consider e.g. the function S2) = {si + Wi, S2 +Pi)- Here for any A > 1, 

if se/ A < < Asi, i = 1,2, then (si -|- Wi)/A < s'l + Wi < A(si + Wi) and 
(s2 +Pi)/A < s'2 + Pi < A{s2 + Pi), therefore Fi{si,S2) is (D,A)-close to 
Fi{s'i,S2). Besides that, adding a constant to si does not change the order 
^qua- The functions id and Hi are treated analogously. 

The other problems considered in Section [2731 either do not admit FPTAS 
unless P=NP (the TSP), or they are solvable in time which is polynomially 
bounded in the input length and thus do not require FPTAS (the SSSP and 
the APSP problems). 
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4.1 Fully polynomial- time approximation scheme 



To identify subsets of states which are {D, A)-close to each other, the al- 
gorithm DPa employs a partition of the set of states into A-boxes (defined 
below). This partition allows to discard "close" states analogously to dis- 
carding of (1 + £:)-dominated solutions which is used in multi-objective opti- 
mization for approximation of Pareto-set (see e.g. 22). The main difference 
is that in our case the states are compared on the basis of their components, 
rather than the components of the vector of objectives. Note that usage of a 
quasi-linear order ^qua in DPa will make {D, A)-closeness only a necessary 
condition for discarding states from consideration. 

Let L be a sufficiently large value, chosen for x and for any required 
precision e G (0,1) (a specific definition of L will be discussed later). To 
describe the algorithm DPa let us consider a family of parallelepipeds that 
constitute a partition of the set B{L, A) = n [0, A^]^: 

{%!,.. .,fc^) ■.h = 0,...,L, £ = 1, . . . , /?}, 

where B(^ki,...,ki3) contains all integer points S = (si, ...,5/3) G Nq, such that: 

( 0, ifA;, = 0, 

se G < [A''^-\ - 1], if < < L, (10) 
\ [A'=^-i,A*^^], ifA;^ = L, 

for all £ G £1 and 

for all ^ G £o- Thus defined parallelepipeds are called A-boxes below. 

Algorithm [3] was suggested in [43j where it was proven to constitute an 
FPTAS with A and L chosen as follows 

A = l + ^, (11) 
2777, 

7ri(n,log2 |x|) ] 

In A ■ ^ ^ 

Equations f lTT]) and f lT2|) ensure L is polynomially bounded in size of the 
input and in 1/e. 



4.2 Fully Polynomial-Time Randomized Approxima- 
tion Scheme 

A family of randomized algorithms over all factors < e < 1 with polynomi- 
ally bounded running times in problem input size |x| and in 1/e that com- 
putes (1 -f £)-approximate solutions with probability at least 3/4 is called a 
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Algorithm 3 DPa for 11 

1: To^So 

2: for i = 1 to n do 

3: 71^0 

4: for S G Ti-i and F G J^j do 

5: let B(^ki,...,kis) be the A-box containing F{S) 

6: if < and ^5' G Tin F{S) ^q,a ^' then 

7: 7- ^ (7- \ G 7- n I 5' F{S)}) U 

8: end if 

9: end for 

10: end for 

11: find S* eTn such that 

r'(Q*\ — / : S G Tn} in case of minimization, 

max{G(S') : 5* G Tn} in case of maximization 

12: return y{S*) 



fully polynomial-time randomized approximation scheme (FPRAS) [23] . The 
constant 3/4 in the definition of FPRAS for optimization problems may be 
replaced by any other constant from the interval (0,1). 

The DP-based EA framework proposed in Section [3] may be modified to 
obtain an evolutionary FPRAS for DP-benevolent problems. 

Now a new relation is defined to substitute ^dom in Algorithm [2j Let 
us introduce the following relation: {i,S) {i',S'), iff Hi{S) > or the 
following three conditions hold: 

a) i = i' 

b) there exist such /ci, . . . , A;^ that S, S' G B[ki,...,kp) 

The EA using this relation is denoted EAa in what follows. 

For an arbitrary S E Si let 6{i,S), i = l,...,n, be the first iteration 
number, when an individual {i,T) was added into population, such that: 

(i) T is {D, A*)-close to S and 

(ii) S ^qua T. 
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In all iterations following 6{i, S) the population will contain an individual T 
that satisfies the conditions (i) and (ii) as well. 

The following lemma indicates that for any non-dominated S & Si, in a 
number of iterations that is on average polynomially bounded in |x| and 1/e, 
an individual {i, T) will be obtained such that T is {D, A*)-close to S and 
•S* ^qua T . The proofs of the lemma and the theorem below are provided 
in [T^ but since this publication might be difficult to access, we reproduce 
the proofs here. 

Lemma 1. Let 11 he DP-benevolent with dimension 13. Then for any stage 
i = 0, . . . ,n, any non- dominated state S in Si and L chosen as defined in 
Equation EH it holds that 

i 

E[e{i, S)] < n{Ln2{n, log, \^\))^ ■ 

k=l 

Proof. Let us use induction on i. For i = the statement holds trivially. 
Consider any state S which is non-dominated in Si. Suppose i > and the 
statement holds for i — 1. 

Lemma 4.7 in [13] implies that there exists a state S** non-dominated 
in Si-i and a mapping F* e J^i, such that F*{S*) = S. Note that the 
induction hypothesis gives an upper bound on expectation of 9{i — 1,5**), 
which is the expected number of iterations until an individual {i — 1,T*) is 
obtained, such that T* is (D, A*~^)-close to 5* and ^qua T*- Again, 
let the mutation that applies F* to an individual {i — 1,T*) be called a 
successful mutation. 

In view of C.2' condition, 

H,{F*{T*)) < H,{S) < 0, 

and by C.l and C.l', either (a) F*{T*) is (D , A'-^)-c\ose to S and S ^q^a 
F#(T#),or (b)5^dom F#(T#). 

In ), after a successful mutation, the population will contain the 

element (i, F*(T*)), or some other element {i, T') such that T' belongs to the 
A-box i3(fci,...,A:^), which also contains F'^iT'^) and besides this F'^iT'^) :<qua. 
T' . After this mutation the population will contain an individual (i, T), such 
that T is (D, A)-close to F*iT*) and F*{T*) ^q^a T. Now since F*iT*) 
is (D, A*~^)-close to S, by the definition of closeness, T is (D, A'')-close to 
S. Besides that, S ^qua F'^iT'^) :<qaa. T, consequently, 5* ^qua T. Thus, in 
case ^a ), successful mutation ensures presence of the required representative 
for S in population on stage i. 
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In case (b), a successful mutation will yield the individual (i, S), since 5* 
is a non-dominated state, and 5* ^dom -F*(T*). After such a mutation, the 
population will contain an individual {i, T), such that T is {D, A)-close to S 
and 5* ^qua T. Obviously, T is also {D, A*)-close to S then. 

To complete the proof it remains to estimate the expected number of 
mutation attempts 9* until a successful mutation occurs, conditioned that 
an individual a* = {i — 1,T*) belongs to the current population V. Note 
that the probability of a successful mutation is 

p* = [n-\{{i-l,S')eV}\-\Hy\ 

at the same time, 

n n 

\V\ = S') eV}\<Y, im, n S, + 0}|. (13) 

i'=\ i'=l 

Consider a single term in the right-hand side of inequality f fT3|) with any 
fixed i' . For each I ^ L\ the index ki may take at most L different values. 
Besides that, in view of condition C4 (iv), for each G £o the coordinate s^i 
characterizing the states from the set Sii may take at most 7r2(n, log2 |x|) 
values. 

Thus, the right-hand side of inequality (|T3l) can not exceed 

nLl^il7r2(n,log2 |x|)l^«l < n(L7r2(n, log2 |x|))'^. 

The statement of the lemma for phase i follows from the fact that E\§{i^ S*)] = 
- 1,T#)] + □ 

The bound on E^{%, 5*)] obtained in Lemma [1] is used to choose the 
stopping criterion for the algorithm EA^. Let the algorithm terminate after 

n 

T = 4n(L7r2(n, log2 |x|))^ ^ \J^,\ (14) 

1=1 

iterations. 

Theorem 3. If the problem U is DP-benevolent, then the family of algorithms 
EAa where A and L are chosen according to / fii]) and KW^. using the stopping 
criterion ([7^ gives an FPRAS. 

Proof. In view of ([H]) and C3 (ii), there exists a non-dominated state 
S* G Sn, such that OPT{^x) = G{S*). By Lemma [H on average within at 
most n(L7r2(n, log2 |x|))^ ■ J2^=i iterations of EAa, a population will be 
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computed, containing an individual {n,T*), such that T* is {D, A'^)-close to 
S* and S* <,ua T*. 

Let us first consider the case where 11 is a minimization problem. By 
condition C3 (i): 



The latter inequality follows from the observations that 7n > 1, (1 + 2^)'^" 
is a convex function in e on the interval e G [0, 2], and the indicated inequal- 
ity holds for both endpoints of this interval. In the case of maximization 
problem 11 analogously we obtain G(T*) > (1 + £:)~^OPT(x). 

Finally, by means of backtracking, a (1 + £)-approximate solution y{T*) 
may be computed efficiently. 

Execution of EA^, with stopping criterion f[T^ . according to the Markov 
inequality, does not yield a (1 + £)-approximate solution with probability at 
most 1/4. 

Finally, by condition C.4, the runtime of each iteration of the EAa is poly- 
nomially bounded in the input length and in 1/e. Summing up the observed 
facts, we conclude that the proposed family of the algorithms constitutes an 



5 Conclusions 

We have examined how to choose a representation for an evolutionary algo- 
rithm such that it obtains the ability to carry out dynamic programming. 
Based on a general framework for dynamic programming we have given a 
framework for evolutionary algorithms that have a dynamic programming 
ability and analyzed the optimization time of such an algorithm depend- 
ing on the corresponding dynamic programming approach. By considering 
well-known combinatorial optimization problems, we have shown that our 
framework captures most of the known DP-based evolutionary algorithms 
and allows to treat other problems. 
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A Bellman Principle for Single-Objective Prob- 
lems 

In this appendix, we describe the Bellman optimality principle in terms of 
the DP method defined by recurrence ([T]). Consider a single-objective maxi- 
mization problem 11. Let 2 < /3 G N be a constant and S C M^, so that the 
first component Si of a state S E S characterizes a quality of the state in 
some sense. 

The Bellman principle applies to a DP algorithm for 11 if the following 
statement holds. Suppose that starting from some state Sq G Sq, a sequence 
of "decisions" Fi E J-'i, . . . , Fn E J-'n leads to an optimal solution for 11. Let us 
denote S* = Fj(Fj_i(...Fi(S'o)...)) €5,, i = 1, . . . ,n. Then for any particular 
state S* = (s^j, . . . , s^j), the subsequence Fi, . . . , Fj is an optimal policy for 
reaching the set of states coinciding with S* in components S2, ■ ■ ■ , sp. By 
an optimal policy here we mean that for any sequence F{ G J-*!, . . . , -F/ G J^i 
starting with some 5*0 G Sq, such that S*^ = F^(F^_;^(...F/(S'o)...)) G Sk, k = 
and S'i = (4, s*2i, s^J, holds s[i < sl^. 

If the Bellman principle applies to a DP algorithm, then for any S2, ■ ■ ■ ,3^ 
it is possible to keep only one state which dominates all states in the subset 
{S' E Si : S2 = S2, . . . , = s^} without a risk to loose optimality of the 
DP algorithm result. 

B Genetic Mechanisms Corresponding to the 
Mutation Proposed in the EA 

The mutation operator proposed in the EA in Section [3^ is a special case of 
the point mutation, where a gene Ai subject to change is selected as the first 
gene which has never been mutated so far (i.e. the first "undefined" gene). 
Such type of mutation may be imagined in a biological system as follows. 

Suppose that for each phase i, i = l,...,n, there is a "controlling" 
gene 5j. The required localization of mutations in gene A^, when Ai is 
the first "undefined" gene, is caused by insertion of some mobile DNA se- 
quence Ci (e.g. a transposon, see [ID]), that can enter the locus of gene Ai, 
and only this locus. We can additionally assume that a mobile element Ci 
is produced if and only if the gene Bi is active (i. e. Bi is subject to tran- 
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scription in the parent individual). Besides that, we can assume that gene Ai 
in the "undefined" condition is silencing the transcription of gene -Bj+i, but 
any mutated state of gene Ai activates the transcription of gene -Bj+i and 
silences the gene Bi. 

Then one can assume that in the i-th generation, i = 1, . . . ,n, only the 
gene Bi is active among Si, . . . , provided that initially only the gene Bi 
was active. At the same time, in the i-th generation, i = l,...,n, the 
insertion mutations occur only in the gene Ai. 

In nature, an example of a mutually exclusive genes activation is observed 
in malaria parasite Plasmodium falciparum. The transitions from one variant 
of a gene to another one depend on the currently active gene variant [23]. 
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